
#include "code.h"

double * polyfit(double * x_input , double * y_input , double * w_input , double power_input , double num_point_input){/*
主体函数，应该包含：
    // 1. (x , y)坐标值
    // 2. 权重
    // 3. 拟合方程次幂
    // 4. 输入坐标的个数
*/

    
    inital_data(x_input , y_input , w_input , power_input , num_point_input);


    return my_algorithm();


}

void mydebug(){
    int i , j;
    printf("\n");
    for ( i= 0 ; i < order ; i ++ ){
        for (j = 0 ; j < order ; j++ ){
            printf("| %le " , G_matrix[i][j]);
        }
        printf ("|     | %le | \n" , d_matrix[i]);
    }
}

double calculate_in_G( int index_r , int index_c){
    int i , j;
    
    int power = index_r + index_c;
    //决定x的次幂数


    if(power == 0){//如果次幂为0 ， 则将权值全部相加即可
        double sum = 0.0;
        for ( i = 0 ; i < num_point ; i++ ){
            sum += weight_per_point[i];
        }
        return sum;
    }
    

    double tmp , sum;
    //tmp为每一个x坐标的次幂数
    //sum为求和


    sum = 0.0;
    for ( i = 0 ; i < num_point ; i++){
         tmp = 1.0;
         for ( j = 0 ; j < power ; j++ ){
             tmp *= x_array[i];
         }//次幂

         sum += weight_per_point[i] * tmp;//求和,要带加权
    }


    return sum;
}

double calculate_in_d( int index_x ){
    int i , j;
    double sum,tmp;


    if ( index_x == 0 ){
        sum = 0.0;
        for ( i = 0 ; i < num_point ; i++){
            sum += weight_per_point[i] * y_array[i];
        }
        return sum;
    }

    sum = 0.0;
    for ( i = 0 ; i < num_point ; i++ ){
        tmp = 1.0;

        for ( j = 0 ; j < index_x ; j++){
            tmp *= x_array[i];
        }

        sum += weight_per_point[i] * tmp * y_array[i];
    }
    
    return sum;

}

int get_max_item_c(int col_num , int s_row_num){//寻找每一列的最大元素
    int o_row_num;
    int i;
    double max_element = -111111.1;
    for ( i = s_row_num ; i < order ; i ++){
        if ( G_matrix[i][col_num]  > max_element ){
            max_element = G_matrix[i][col_num];
            o_row_num = i;
        }
            
    }

    return o_row_num;
}

void swap_matrix_line(int line_2_swap1 , int line_2_swap2){//交换指定两行
    int i = line_2_swap1 , j = line_2_swap2;
    int k;
    for (k = 0 ; k < order ; k++){
        double tmp = G_matrix[i][k];
        G_matrix[i][k] = G_matrix[j][k];
        G_matrix[j][k] = tmp;
    }
    
    double tmp2 = d_matrix[i];
    d_matrix[i] = d_matrix[j];
    d_matrix[j] = tmp2;


    // mydebug();
    //check

}

double my_round(double num){//用于将过小的数字约为0
    double temp = num;
    if (temp < 0)
        temp = -temp;
    
    
    return temp < 1e-7 ? 0.0 : num;//如果数字小于 1e-7，则认为由于精度问题，导致浮点数无法靠近0
}

void calculate_row_x(int row , int col){//计算第x列的主元素分析
    int i , j;

    double tmp = 1.0 / G_matrix[row][col]; // 因子，用于消元

    for ( i = row+1 ; i < order ; i++ ){
        double tmp2 = G_matrix[i][col] * tmp;
        d_matrix[i] = my_round( d_matrix[i] - d_matrix[row] * tmp2 );
        for ( j = col ; j < order ; j++ ){
            G_matrix[i][j] = my_round (G_matrix[i][j] - G_matrix[row][j] * tmp2);
        }
    }
    // printf("---------in----------");
    // mydebug();


}

double * resolve_matrix(){
    int i , j;
    double * answer;
    answer = (double*)malloc(order);
    for (i = 0 ; i < order ; i++){
        answer[i] = 0;
    }

    for ( i = order - 1 ; i >= 0 ; i--){
        double tmp = d_matrix[i];
        for (j = i ; j < order ; j++){
            tmp = tmp - answer[j] * G_matrix[i][j];
        }
        answer[i] = my_round(tmp / G_matrix[i][i]);
        printf("| ---answer%d = %le--- " , i,answer[i]);
    }
    printf("\n");
    return answer;
}

void calculate_matrix(){
    int i , j ;
    for ( i = 0 ; i < order ; i++ ){
        int max_line_num = get_max_item_c(i , i);
        swap_matrix_line(i , max_line_num);
        mydebug();
        calculate_row_x(i , i);
    }
    // mydebug();



}

double * my_algorithm(){
    

    int i , j;
    for ( i = 0 ; i < order ; i++ ){
        for ( j = 0 ; j < order ; j++ ){
            G_matrix[i][j] = calculate_in_G(i , j);
            printf("| %le " , G_matrix[i][j]);
        }
        d_matrix[i] = calculate_in_d(i);
        printf ("|     | %le | \n" , d_matrix[i]);
    }//计算两个矩阵
    //check

    calculate_matrix();

    return resolve_matrix();

}

void inital_data(double * x_input , double * y_input , double * w_input , double power_input , double num_point_input){
    // 输入应该包括：
    // 1. (x , y)坐标值
    // 2. 权重
    // 3. 拟合方程次幂
    // 4. 输入坐标的个数

    int i;

    order = power_input + 1;
    num_point = num_point_input;

    for ( i = 0 ; i < num_point_input ; i++){
        // printf("flag\n");
        x_array[i] = x_input[i];
        y_array[i] = y_input[i];
        weight_per_point[i] = w_input[i];
    }
    
}
